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The ionization of a one-dimensional model Helium atom in short laser pulses using time-dependent 
density functional theory is investigated. We calculate ionization probabilities as a function of laser 
intensity by approximating the correlation function of the system adiabatically with an explicit 
dependence on the fractional number of bound electrons. For the correlation potential we take the 
derivative discontinuity at integer numbers of bound electrons explicitly into account. This approach 
reproduces ionization probabilities from the solution of the time-dependent Schrodinger equation, 
in particular the so-called "knee" due to nonsequential ionization. 
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Experimental double ionization yields for Helium 
atoms in laser fields of A = 614nm 0] and A = 780nm 0] 
show a "knee"- structure as a function of intensity, i.e., 
double ionization is enhanced over the expected yields 
for sequential ionization by several orders of magnitude 
and was therefore named nonsequential double ionization 
(NSDI). 

Theoretically, the problem was addressed successfully 
using the strong-field approximation (see, e.g., 0] and 
references therein) and classical methods |4| . The widely 
accepted mechanism behind NSDI relics on the rescatter- 
ing of the first electron with its parent ion, collisionally 
ionizing (or exciting) the second electron. 

However, investigations using time-dependent density 
functional theory (TDDFT) [jf failed to reproduce the 
characteristic "knee" of the double ionization probability 
0,0,0 ■ This failure is especially important given the fact 
that TDDFT offers the possibility to treat the strong- 
field dynamics of atomic and molecular systems from first 
principles in a numerically tractable way, whereas the so- 
lution of the Schrodinger equation in full dimensionality 
0] is limited to two active electrons. Other methods used 
to tackle the NSDI problem from first principles are R- 
Matrix Floquet theory |lC Lor multi-configurational time- 
dependent Hartree-Fock n], which, however, have 
not yet reproduced NSDI-" knees" either. 

Two reasons for the failure of TDDFT in the context 
of NSDI were identified ^jjj : The approximation of the 
exchange-correlation functional u xc and the approxima- 
tion of the functionals for the calculation of the ionization 
probabilities. Initial attempts to overcome these prob- 
lems were not successful [l3|. Recent progress in ap- 
proximating v xc has been made by including the deriva- 
tive discontinuity at integer numbers of bound electrons 
[Til Hil . In this work we therefore focus our atten- 
tion on the open problem of constructing functionals for 
the calculation of the ionization probabilities. We show 
that by using a rather simple adiabatic approximation of 
the correlation function, the ionization yields for a model 
He-atom in few-cycle laser pulses (for which we have the 



exact time-dependent Schrodinger solutions at hand) are 
remarkably well reproduced. 

We investigate Helium in linearly polarized N — 3-cycle 
laser pulses of duration T = N 2 tt/u; with a sin 2 -pulse en- 
velope in dipole approximation, i.e., the vector potential 
is of the form A (t) = Asm 2 (^t) sin(wi) for 0<t<T. 
Two differentpulses with frequencies used in experimen- 
tal work 0, 0] are chosen: u = 0.058 and w = 0.074 
(atomic units are used throughout this work unless oth- 
erwise indicated). The linear polarization of the laser 
pulses allows to model Helium by a one-dimensional atom 
with soft-core potentials for the Coulomb interactions. It 
is known that the essential features of the nonsequential 
double ionization process are described well by this model 
EE BUI 0, Hill- For different effective peak inten- 
sities I = 1(A) of the laser pulses we solve the time- 
dependent Schrodinger equation (TDSE) \d t ip = Hip 
with the Hamilton operator 



H=J2 {-\dl k +V{x k ,t)\ + W(x 1 ~x 2 ) , (1) 

fe=l,2 ^ ' 

where H = H (xi, X2,t), the external potential V (x, t) = 
iA{t)d x — 2/Vx 2 + e C c (the A 2 -term has been trans- 
formed away) and the electron-electron interaction po- 
tential W (x) = \j \Jx 2 + e cc . The solution is the two- 
electron wavefunction ip (xi, X2, t) in a spin singlet state. 
Therefore we have just one Kohn-Sham orbital fulfill- 
ing \dt4> — Hks4> with (f> — tfi(x,t) and Hjts(x,t) = 
— (1/2) d 2 + V (x, t) + Uhxc (x, t). The Hartree-exchangc 
part tJhx = + v-x. of the potential «hxc = Vh + + 
v c is known exactly for the Helium atom, t>h (x, t) — 
J dx' 2 \4> (x', t) | 2 W (x - x') and v x (x, t) = ~\ v h (x, t). 
A simple model for the correlation potential which takes 
into account the derivative discontinuity at integer num- 
bers of bound electrons is the recently proposed (l^j 



,,LK05 
J hxc 



[x,t) = 



N /N(t) 



1 + exp [C(N /N(t)-2)] 



Vhx(x,t) , (2) 



where N(t) = dx 2 \(j> (x, t) | 2 is the number of bound 



2 



electrons and No = N(t = 0). C is a positive, sufficiently 
large number, e.g., C = 50. We choose e cc = 0.5 which 
yields the correct 1^ =2.0. In order to reproduce Ip — 
0.904 in the Schrodinger equation e cc = 0.329 and in the 
TDDFT calculation e cc = 0.343 are used (results do not 
change qualitatively for a different choice of e ce and e ee ). 
To encompass the bound states, the parameter a is set 
to a — 6 throughout this work. The Schrodinger and 
Kohn-Sham equations are solved numerically by a split- 
operator time propagator on a numerical grid (see, e.g., 
|l8| and references therein). 

For our model atom the diagonal of the two-electron 
density matrix is the pair density p(x\,X2,t) = 
2 \if> (xi, X2, t) | 2 . The exact density of the system 



n(x,t)= / dx'p(x,x',t) 



(3) 



can in principle be calculated from the Kohn-Sham or- 
bital of a TDDFT calculation with the correct exchange- 
correlation potential u xc as n(x,t) = 2|0(x,i)| 2 . Us- 
ing v xc in the LK05-approximation yields approxi- 
mate densities «lko5- The time-dependent correla- 
tion function of the model atom is g xc [n] (xi , x-x , t) — 
p (xi, X2, t) jn (xi, t) n (x2, t). The conditional probabil- 
ity to find at time t an electron at X2 in dx2 if there is 
an electron at x\ is given by n(x 2 ,t) g xc (xi,x 2 ,t) dx2 



e_p 

5] 



theorem [5j assures that g xc can in principle be expressed 



Since the probability is an observable, the Runge-Gross 

ncipL 

exactly as a functional of the density |l3| . 

For a system with a two-electron wavefunction that can 
be expressed as a product of single-electron orbitals, g xc 
simplifies to g xc — g x — \- The correlation contribution to 
the correlation function is thus defined as 



g c [n] (xi,x 2 ,i) 



p{xi,x 2 ,t) 1 
n (xi, t) n (x 2 ,t) 2 



(4) 



Due to Q, g c fulfills the integration constraints 
/ dxin(xi,t) g c [n] (xi,Xj-n,t) = with i,j e {1,2}. 

A quantity of central importance to our work is the 
number of electrons bound to the Helium core at time t, 



N(t) 



+a 



dx n (x, t) 



(5) 



which also plays a central role in the exchange-correlation 
potential u xc in ionization processes, as in Eq. J5J. 

In the two-electron space of the model atom we at- 
tribute areas to single ionization P + (either |xi| >a and 
|a?2 1 _ a or |xi| < a and \x 2 \ > a) and double ionization 
P 2+ (|xi| > a and |x2| > a). Integration of \ip (xi, X2, T) | 2 
over these areas yields the ionization probabilities. This 
procedure to calculate ionization probabilities is well- 
established 0, 0, and leads to ionization probabilities 
depicted in Fig.^ Substituting the density and the cor- 
relation function for the pair density in the integrals and 



using P" +P+ + P 2+ = 1 (P° being the probability for 
the atom to remain neutral) yields the exact ionization 
probabilities 0] 



P+(t) = 2p(t)[l-p(t)]-I c (t) 
P 2+ (t) = [1- P (t)] 2 + I/ C (t) , 



(6) 
(7) 



where p(t) — ^ J*_^° dx n (x, t) and the correlation inte- 
gral is given by 



Ic(t) 



+a 



dxi / dx 2 n(xi) n(x 2 ) 3c (xi,x 2 ) , (8) 



with the density and the correlation function depending 
on time t and g c — g c [n] . 

It is important to note that Eqs. ©-© assure that 
only the exact density in the range —a < x < +a and 
the exact correlation function for —a < x\,x 2 < +a are 
required to calculate the exact ionization probabilities. 

The functional dependence of g c on n is not known 
and the calculation of the pair density p at time t = T is 
computationally demanding even for the one-dimensional 
Helium model atom. Therefore, it is necessary to find 
a suitable approximation for g c [n](x\,x 2 ,T) to calcu- 
late I c . Setting I c (t) — corresponds to assuming the 
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FIG. 1: Single and double ionization probabilities for exact 
and LK05 densities using J c = compared to the TDSE so- 
lution for a A = 780 nm laser pulse. 

two-electron wavefunction to be a product of identical 
single-electron orbitals. This interpretation of the Kohn- 
Sham orbitals as electron orbitals has been widely used 
0, EH Q|- However, this approach leads to ioniza- 
tion probabilities which differ significantly from the exact 
probabilities even when the exact density (J3J) is used 0] , 
as shown in Fig.^ The LK05 approximation for v xc leads 
to probabilities which closely resemble the probabilities 
calculated from the exact densities. 

The TDSE solution yields the exact pair density. This 
allows us to calculate the exact g c and I c as a reference 
for approximations of g c . Figure reveals a structure of 
the integrand of J c (T) which is remarkably simple. This 
is caused by the multiplication of g c by the density at 
Xi and at x 2l which exhibits a maximum at the origin. 
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FIG. 2: Comparison of the integrand of the exact 7 C (T) (left) 
and the adiabatic approximation 7* (T) (right) for different 
effective peak intensities of a A = 780 nm laser pulse. 



For the other laser pulse the integrand of I c reveals a 
qualitatively similar structure. 

The values of I c (T) for the two laser pulses are de- 
picted in the left panel of Fig.[3]as a function of the num- 
ber of bound electrons N(T). I c shows a quantitatively 
similar behavior for both laser pulses: a minimum value 
for I c is reached when the number of electrons is equal to 
one. This result points to the importance of the quantity 
N(T) for the value of I c (T). 

We shall now introduce our adiabatic approximation 
for g c . The starting point is our assumption that the sys- 
tem in the vicinity of the Helium core after the laser pulse 
closely resembles a groundstate system with a fractional 
number of electrons given by the number of bound elec- 
trons N (T). We approximate numerator and denomina- 
tor of the correlation function separately as linear combi- 
nations of the groundstate densities rii and groundstate 
pair densities pi of i-electron systems. The density of the 
model Helium atom is expressed as 01 



(1 - N) n + Nni 

(2 - N) m + (N- 1) n 2 



<N< 1 

1 <A< 2 



(9) 



where n 



n A (x,t), rii — rii(x) and N — N (t) is 



given by JSJ. Using n A (t) in the groundstate Kohn-Sham 
equation to calculate v c (t) and consequently employing 
v c (t) in the time-propagation has been shown to improve 
ionization probabilities compared to v c = [lfij . Since 
n (x) dx is the probability to find an electron at x in dec, 
the density of the zero-electron system vanishes: uq = 0. 
Assuming the same adiabatic dependence of p A (t) on 



P A = 



(1 

(2 



N) po 
N) Pl 



Npi 
(N- 



1 



Pi 



<N< 1 

1 <7V< 2 



(10) 



where p A — p A (x\, X2,t), pi — pi {x\, x 2 ) and N = N (t). 
The pair density p{x\,X2) gives the joint probability to 
find an electron at x\ in dx\ and an electron at x 2 in 
dx 2 ■ Corresponding to our consideration for the density, 
the pair densities of the zero- and one-electron systems 
are zero: po = 0, p% = 0. For our model atom, m (x) 
is the groundstate density of He + . It is calculated from 
the solution of the respective one-dimensional stationary 
Schrodinger equation with e cc = 0.5. Consistently we use 
n 2 (%) = n (x, t = 0) where n denotes the density utilized 
to calculate N(t) in (|SJ). Our adiabatic approximation of 
g c thus is 



1 
1 



p A (xi,X2,t) 

(xi,t) n A (x 2l t) 



1 
5 



< iV< 1 
1< A^< 2 



(11) 



with g A = g A (xi,x 2 ,t). Although for N(T) = 1 one has 
g A = 1/2 which yields the unphysical value P 2+ (T) = 0, 
good results are already acquired for N (T) only slightly 
smaller or larger than one, as is shown below. 

In this form the approximation recovers the exact cor- 
relation function for the groundstate. This approxima- 
tion can only be valid in the immediate vicinity of the 
Helium core. However, as pointed out above this is suffi- 
cient to calculate the correlation integral I A (T) . The in- 
tegration constraints J dxi n (xi, t) g c [n] (xi, x^ ut) = 
with i,j £ {1,2} arc not fulfilled since for \xi\ > a we 
have n/0 and g c ^ but set g A = 0. The same densi- 
ties as in Eq. © are used in the integrand of I A . 
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FIG. 3: Value of the exact I c (T) as a function of the number 
of bound electrons (left). Value of 7^ (T) using exact densities 
compared to the exact 7 C (T) (right). Results for both laser 
pulses are shown together with differently shaped symbols. 

The groundstate inputs needed to calculate g A are ac- 
cessible numerically also for more complex systems than 
Helium. Moreover, in the three-dimensional case the 
groundstate pair density can be calculated approximately 
via one of the well-known approximations for the ground- 
state correlation function |19l |2fJ . 
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The right hand side of Fig.0 shows the integrand of 
1^ (T) using the exact densities. The qualitative agree- 
ment with the exact integrand of I c (T) is good. Fur- 
thermore, the maximum and minimum values of the in- 
tegrand are close to the exact values. This is a strong 
indication that the adiabatic approximation of g c cap- 
tures the essential features which contribute to I c . In 
the right plot of Fig.|2| the adiabatic approximation l£ 
is compared to the exact values I c for both laser pulses 
as a function of N(T). The adiabatic approximation ex- 
hibits a very similar dependence on N as the exact values. 
However, for N (T) w 1.0 the values differ from the exact 
results. In our adiabatic approximation for N{t) < 1 all 
correlation effects vanish, as can be seen from (|11[) . This 
might indicate the necessity to include memory effects in 
a refined approximation of p. 
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FIG. 4: Single and double ionization probabilities for exact 
and LK05 densities using the adiabatic approximation of the 
correlation integral 1^ compared to the TDSE solution for a 
A = 780nm laser pulse. 

Figure 0| shows the ionization probabilities calculated 
from the exact and the LK05 densities using the approx- 
imated i~ (T) for the A = 780 nm laser pulse. Single 
ionization probabilities agree very well with the exact 
probabilities from the TDSE for both exact and LK05 
densities. The systematically too low values for I c = 
are corrected. Double ionization probabilities show a 
good agreement for high intensities. For low intensi- 
ties the double ionization probabilities (and 1^) are so 
small that already small deviations of 1^ from the ex- 
act I c lead to relatively large deviations in the logarith- 
mic plot of the ionization probabilities and to minima 
and maxima in the double ionization probability. Since 
<7c = 1/2 at N (T) = l, the adiabatically corrected prob- 
abilities for the two laser pulses examined show a min- 
imum at intensities when the second derivative of the 
exact probabilities changes its sign at the onset of the 
"knee" -structure (where N (T) ~ 1). Comparison with 
the values for J c = (Fig.QJ evidences a significant im- 
provement of the TDDFT ionization probabilities using 
our approximation of the correlation function l|ll(l . 

In summary, to acquire the exact, multiple ioniza- 
tion probabilities of atoms in strong laser fields from 
time-dependent density functional calculations, knowl- 



edge of the correlation integral I c (T) is vital. For a 
one-dimensional Helium model atom, an adiabatic ap- 
proximation of the correlation function g c close to the 
Helium core yields an approximation of I c which agrees 
qualitatively and quantitatively well with the exact I c . 
The exact single ionization probabilities and the "knee"- 
structure in the double ionization probability are repro- 
duced. The (pair) density of Helium after the laser pulse 
can be approximated by a linear combination of atomic 
and ionic groundstate (pair) densities using the fractional 
number of bound electrons. This method opens the pos- 
sibility to apply the adiabatic approximation in three di- 
mensions and to the ionization of more complex atoms. 

We would like to thank Manfred Lein for stimulating 
discussions and the Deutsche Forschungsgemeinschaft for 
supporting this work. 
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